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ABSTRACT 

We have numerically explored the stable planetary geometry for the multiple 
systems involved in a 2:1 mean motion resonance, and herein we mainly concen- 
trate on the study of the HD 82943 system by employing two sets of the orbital 
parameters (Mayor et al. 2004). In the simulations, we find that all stable orbits 
are related to the 2:1 commensurability that can help to remain the semi- major 
axes for two companions almost unaltered over the secular evolution for 10 7 yr, 
and the apsidal phase-locking between two orbits can further enhance the sta- 
bility for this system, because the eccentricities are simultaneously preserved to 
restrain the planets from frequent close encounters. For HD 82943, there ex- 
ist three possible stable configurations: (1) Type I, only 9\ « 0°, (2) Type II, 
0! « 9 2 « #3 « 0° (aligned case), and (3) Type III, 6 X w 180°, 9 2 « 0°, 9 3 « 180° 
(antialigned case), here the lowest eccentricity- type mean motion resonant ar- 
guments are 6i — X% — 2A 2 + W\ and 9 2 = X\ — 2A 2 + vc 2 , the relative apsidal 
longitudes 9 3 = w\ — vj 2 = Aw (where Xi j2 are, respectively, the mean longitudes 
of the inner and outer planets are the longitudes of periapse). And we find 
that the other 2:1 resonant systems (e.g., GJ 876 or HD 160691) may possess one 
of three stable orbits in their realistic motions. In addition, we also propose a 
semi-analytical model to study ej — Aw Hamiltonian contours, which are fairly 
consistent with direct numerical integrations. With the updated fit, we then ex- 
amine the dependence of the stability of this system on the relative inclination, 
the planetary mass ratios, the eccentricities and other orbital parameters: in the 
non-coplanar cases, we find that stability requires the relative inclination being 
~ 25° or less; as to the planetary mass ratio, the stable orbits for HD 82943 
requires sini > 0.50 for a fixed value or mi/m 2 < 2 where m 2 remains for the 
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varying mass ratio; concerning the eccentricities, the system can be always steady 
when < e<i < 0.24 and < e\ < 0.60. Moreover, we numerically show that the 
assumed terrestrial bodies cannot survive near the habitable zones for HD 82943 
due to the strong perturbations induced by two resonant companions, but these 
low-mass planets can be dynamically habitable in the GJ 876 system at ~ 1 AU 
in the numerical surveys. Finally, we present a brief discussion on the origin of 
the 2:1 resonance for HD 82943. 

Subject headings: methods:N-body simulations — celestial mechanics — plane- 
tary systems — starsdndividual (HD 82943, GJ 876, HD 160691) 

1. Introduction 

The discovery of the extrasolar planets is opening a new world beyond our solar system. 
Ever since 1995, Mayor & Queloz (1995) detected the first extrasolar giant Jupiter-51 Peg, 
and to date there are more than 100 planetary systems discovered 1 (Marcy, Cochran, & 
Mayor 2000; Butler et al. 2003) using the radial velocity technique in the surveys of nearby 
young stars. At the time of writing, a dozen of multiple planet systems-HD 82943, GJ 876, 
HD 168443, HD 74156, 47 Uma, HD 37124, HD 38529, HD 12661, 55 Cnc, v And (Fischer 
et al. 2003), HD 169830 (Mayor et al. 2004) and HD 160691 (Jones et al. 2002) were 
discovered in recent years and this number is undoubtedly cumulative as the measurements 
are being carried on. Hence, it is necessary to categorize the discovered multiple planetary 
systems according to their statistical characteristics (such as the distribution of the planetary 
masses, semi-major axes, eccentricities and metallicity) (Marcy et al. 2003), then to study 
the correlation between mass ratio and period ratio (Mazeh & Zucker 2003) and further to 
improve the understanding of the relationship between the planet occurrence rate and stellar 
metallicity (Santos et al. 2003; Fischer, Valenti, & Marcy 2004). The other key point is to 
investigate possible stable configurations for the multiple systems, in which the observations 
reveal that most of them are typically characterized by mean motion resonance (MMR) and 
(or) apsidal phase-locking between their orbiting companions (Fischer et al. 2003; Ji et al. 
2003a; Lee & Peale 2003), so that one can better understand the full dynamics of these 
systems. In the present study, we mainly focus our attention on the HD 82943 system. 

As of April 4, 2001, the Geneva Extrasolar Planet Search Team 2 announced the discov- 



1 see also http://cfa-www.harvard.edu/plancts/bibli.html and http://exoplanets.org/ 
2 see http://obswww.unige.ch/~udry/planet/hd82943syst.html 
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ery of the HD 82943 planetary system at an ESO press release, consisting of two Jupiter-like 
planets orbiting the parent star. Subsequently, Israelian et al. (2001, 2003) reported the 
discovery of 6 Li in the atmosphere of this metal-rich solar-type star and indicated that the 
presence of 6 Li can probably be shown as evidence for a planet or planets having been en- 
gulfed by the star HD82943. At first, let us review several stellar features for this star:HD 
82943 is a GO star with B — V = 0.623, Hipparcos parallax of 36.42 mas, distance of 27.46 
pc, and [Fe/H] = 0.29 ± 0.02 (Santos et al. 2001, 2003). The star is aged 2.9 Gyr with 
the stellar mass of 1.15 M Q . Similar to the GJ 876 system (Marcy et al. 2001), which is 
believed that stability of the planetary system might be sustained by the 2:1 mean- motion 
resonance (MMR) and apsidal alignment in the semi-major axes (Kinoshita & Nakai 2001; 
Lee & Peale 2002; Ji, Li, & Liu 2002), the two planets of HD 82943 are now also close to a 2:1 
commensurability, with orbital periods 435.1 ± 1.4 and 219.4 ± 0.2 d, and semi-major axes 
1.18 and 0.75AU (Mayor et al. 2004). However, the previous investigations (Gozdziewski 
& Maciejewski 2001) of the stability of the HD 82943 system showed that the system given 
by the earlier best-fit orbital solution was extremely unstable due to the strong interaction 
of two massive planets with high eccentricities, hence a natural question is that whether 
the two companions for this system are really in a 2:1 orbital resonance and whether the 
system is secularly stable with the observed best-fit orbital parameters. Additionally, the 
stability for HD 82943 places constraints on both planetary masses-if the mass of the outer 
companion is smaller than that of the inner planet (Hadjidemetriou 2002), a system that 
harbors two planets moving on elliptic 2:1 resonant periodic orbits will be destabilized as a 
finale. Nevertheless, the recently updated orbital solutions (Mayor et al. 2004), with more 
abundant observations as of October 13, 2003, show that the inner planet is a bit more 
massive than the outer one, thus it is indeed noteworthy to thoroughly reanalyze the ob- 
served HD 82943 system, still these mentioned issues of this fascinating system are to arouse 
our great interests and inspire us to seek after the answers. Consequently, our first goal 
is to carefully explore the system to fully understand in what kind of likely configurations 
that two companions may orbit about their host star and if exist, to reveal the dynamical 
mechanisms to remain the system. 

In a viewpoint of celestial mechanics, the former studies are concerned for the dynamical 
analysis for the global stability of the exosystems (Gozdziewski & Maciejewski 2001, 2003; 
Ji et al. 2002; Kiseleva-Eggleton et al. 2002 ; Barnes & Quinn 2004), the investigations of 
resonant picture in the phase space (Hadjidemetriou 2002; Ji et al. 2003a; Haghighipour 
et al. 2003; Callegari Jr., Michtchenko, & Ferraz-Mello 2004) or chaotic behavior (Jiang 
& Yeh 2004) in these systems, and the existence of Habitable Zones (HZ) for Earth-like 
planets that supply stable liquid-water environment to cultivate extraterrestrial intelligent 
beings (Kasting, Whitmore, & Reynolds 1993; Jones & Sleep 2002; Laughlin, Chambers, & 
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Fischer 2002; Dvorak et al. 2003; Menou & Tabachnik 2003). Other recent works on the 
2:1 resonant geometry of the extrasolar planetary systems include those by Hadjidemetriou 
& Psychoyos (2003), Beauge, Ferraz-Mello, & Michtchenko (2003) and Lee (2004). They 
studied the presence and location of stable equilibrium solutions or dynamical evolution of 
such exosystems in more general way. In this work, we revisit the dynamics of the HD 
82943 system and outline several possible stable configurations for a system harboring two 
companions in a 2:1 resonance, and our new findings are that this system can be locked in 
aligned or antialigned orbital motions, which these results are never presented nor analyzed 
in other studies. Specifically, at first, we introduce the numerical setup (see §2) for the 
dynamical simulations. In §3, we make a study of the various planetary configurations for 
HD 82943 by means of the semi-analytical treatment and direct integrations over long-term 
orbital evolution. We present the numerical results of HD 82943 by employing two series 
of the best-fit orbital solutions, and we discover three kinds of stable planetary geometry 
for this system that all are linked to a 2:1 resonance. Furthermore, we advance a semi- 
analytical model that can avoid the difficulties in the perturbation expansions for the larger 
eccentricities and still help to explain the numerical evolutions. Still, we investigate how 
the stability for the resonant topology depends on the planetary mass ratio and the orbital 
parameters on the basis of the new fit (Mayor et al. 2004). In §4, we explore whether there 
exist an Earth-like planet surviving about the Habitable zones for the authentic systems of 
HD 82943 and GJ 876. 

The other extraordinary phenomenon in the exosystems is that many of them can host 
giant Jupiter-like planets with strikingly larger eccentricities, and in Figure 1 is shown the 
distribution of the eccentricities of 111 planets 3 discovered at present day. The figure exhibits 
that more than 50% of the planets have the eccentricities larger than 0.30, and HD 80606 b 
(Naef et al. 2001) can occupy the eccentricity amounting up to 0.93, which indicates that 
such unusual discoveries are quite different from the circumstances in our solar system, where 
most of the major planets revolve around Sun on the near-circular orbits. In addition, from 
the view of the observations, the larger eccentricities seem to favor the planet detectability, 
where the semi-amplitude of wobble velocity K oc M p sini/y / a(l — e 2 ) (with M p <C M c ), 
herein M c , M p , a, e and i are, respectively, the stellar mass, the planetary mass, the orbital 
semi-major axis, the eccentricity and the inclination of the orbit relative to the sky plane, 
requires more massive planetary mass (~ the order of Jupiter's mass), a lower a (e.g., close- 
in giant planets) and a higher e (see Fig. 1). Moreover, many previous authors have put 
forward diverse theories and possible mechanisms to interpret the evolution and origin of the 
orbital eccentricity variation in the exosystems: the disk-planet or planet-planet interaction 



3 The data were taken from http://cxoplanets.org, as of Aug. 5, 2003 
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(Goldreich & Tremaine 1980; Weidenschilling & Marzari 1996; Lin & Ida 1997; Ford, Rasio, & 
Sills 1999) render the orbital migration of gas giant planets embedded in the protoplanetary 
systems (Lin, Bodenheimer, & Richardson 1996; Ward 1997; Bryden et al. 2000; Nelson 
et al. 2000) and the eccentricities are usually believed to excite to the observed values 
through complicated processes (Snellgrove, Papaloizou, & Nelson 2001; Kley 2003) or hybrid 
mechanisms (Lee & Peale 2002; Chiang, Fischer, & Thommes 2002) related to both the 
planetary migration and the resonant capture. Therefore, as a second objective, we intend 
to preliminarily make a discussion of the possible origin of the eccentricities or resonant 
configurations (see §5) for two planets of HD 82943. As a final part, in §6, we summarize 
our principal results and give a brief discussion. 



2. Numerical setup 

In the present work, we aim to numerically investigate the orbital motions for the two 
companions for HD 82943 in a three-dimensional space. And at first, let us bear in mind 
that the two companions of this system are assumed to be in the same orbital plane for our 
simulations except where noted. Here we adopt the best-fit orbital parameters from the web 
site of the Geneva Team. As is known, the two-Kepler fit produce five parameters set-(Ki, 
Pi, ej, Ui, T pi ), where for each planet i — 1,2, are the amplitude Ki, the orbital period Pi, the 
eccentricity e^, the argument of periapse u>i and the time of periapse passage T pi . The Geneva 
Team presented two sets of the orbital data and hereafter we respectively call them Fit 1 and 
Fit 2. Here for Fit 1 of 107 observations (as of July 31, 2002) with the residual of 7.4 m/s, 
the data are listed here: mi sini = 0.88 M Jup , m 2 sini = 1.63 M Jup , a\ = 0.73 AU, a 2 = 1.16 
AU, ei = 0.54 ± 0.05, e 2 = 0.41 ± 0.08, uj ± = 138 ± 13° and oo 2 = 96 ± 7°, where hereafter the 
subscripts 1 and 2 denote the inner and outer planets, respectively. And Fit 2 were derived by 
fitting 142 observations (as of October 13, 2003, see also Mayor et al. 2004) with the residual 
of 6.8 m/s, the orbital parameters are presented: mi sin? = 1.85 Mj up , m 2 sini = 1.84 Mj up , 
ai = 0.75 AU, a 2 = 1.18 AU, e x = 0.38 ± 0.01, e 2 = 0.18 ± 0.04, uji = 124 ± 3° and 
uj 2 = 237 ± 13°. Apparently, the eccentricities, semi-major axes, periapse arguments, true 
anomalies and others can act as key factors to determine the shape of the radial velocity 
curves. However, due to the large interaction between two planets in a low mean motion 
resonance, Laughlin & Chambers (2001) pointed out that short-term perturbations among 
massive planets in multiple planet systems (such as GJ 876, HD 82943) can result in radial 
velocity variations of the central star that differ substantially from velocity variations derived 
assuming the planets are executing independent Keplerian motions. On the other hand, the 
best-fit orbital solutions can vary a little from time to time provided that more updated 
observations were supplied or other elaborate fitting procedure was considered (Mayor et 
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al. 2004), or the velocity trend is removed due to the unveiling of an additional planet 
or more planets (Fischer et al. 2001) for a specific system. Therefore, at this stage, it is 
meaningful to study the dynamical characteristics of the system based on the subsequent 
points: (i) slightly varying the orbital parameters about the best-fit parameters by adding 
the uncertainties with direct integrations, (ii) seeking for the possible stable geometry for HD 
82943 in parameter space, (iii) the dynamical mechanisms of retaining the system and (iv) 
the likely eccentricity or configuration origin for two planets, which motivate us to dedicate 
to this contribution. 

Next, we continue to describe the means of preparation for the initial orbital data before 
we start the numerical integrations: for each planet, we usually need to obtain six orbital 
elements - the semi-major axis a, the eccentricity e, the orbital inclination /, the nodal 
longitude Q, the argument of periastron uj and the mean anomaly M. Here, two groups of 
the initial orbital elements were produced on the basis of the data by Fit 1 and Fit 2. In each 
group, we assume that the semi-major axes of the two planets are always unchanged, say, 
for Fit 2, being 0.75 and 1.18 AU respectively. As we have indicated above, the two resonant 
planets are supposed to be coplanar for the initial configuration and the inclinations are both 
taken as small constant values not far from zero, e.g., 0.5°. The eccentricities and arguments 
of periapse are generated in the orbital parameters space in the proximity of the best-fit 
orbital solutions given the nominal observation errors. For example, for Fit 2, we take the 
observed eccentricities e\ or e2 to be centered, respectively, at 0.38 and 0.18, and randomly 
displaced by not more than 3a, respectively, then we obtain the resulting initial eccentricities 
for integration. Still, we carried out similar steps to achieve the starting values of periapse 
arguments. The remained two angles of nodal longitudes and mean anomalies 4 are randomly 
distributed between 0° and 360°. As a consequence, we obtained 100 orbits for each planet, 
to perform the integration of the planetary system, and each pair of the orbits was integrated 
for the time span of 10 7 yr, unless otherwise stated. Obviously, these orbits just lie in the 
neighborhood of the fitting parameters and represent the orbital motions quite close to the 
reality. Here we make study of the possible motion near the best-fit orbital solutions and 
expect that such investigations can be helpful to reveal some important dynamical features 
for the studied system or present valuable clues on the searching for other exosystems related 
to a 2:1 resonance. 



4 Laughlin G. and Lee M.H. (2003, private communications) pointed out to us that the mean anomalies 
and periapse arguments should perfectly match the radial velocities. Here, we choose proper values for M, 
u and fi, such that the linear combination of initial arguments can be close to the resonant configuration, 
this is reasonable, because the real system is observably stable and approximate to 2:1 MMR. And in Table 
1 and 2, we present the orbital solutions where the measurement errors of Aw = and Ae = 0, indicating 
that the eccentricities and periapse arguments are just equal to those by the Geneva Team. 
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With the N-body codes (Ji et al. 2002), we carried out the numerical integrations for 
the HD 82943 system. In our simulations, we adopted the mass of the host star to be 1.05 
M & for Fit 1, and 1.15 M & for Fit 2. Under the assumption of sim = 1, the minimum 
masses of two companions are employed throughout the paper unless we state elsewhere. 
In addition, we utilized the time step to be one percent of the orbital period of the inner 
planet (HD 82943 c) in the integrations. The numerical errors were effectively controlled 
during the integration, with the local truncation error 10~ 14 for the time span of 10 7 yr, in 
the meantime the accuracy was also checked by the energy errors. In comparison, we also 
employed symplectic integrators (Feng 1986; Wisdom & Holman 1991) to integrate the same 
orbit to assure the results for some cases. 

In the research of the planetary systems, the maximum Lyapunov characteristic expo- 
nent (LCE) is usually adopted to identity the regular or chaotic orbits, since the chaotic 
orbits are sensitively dependent on the initial conditions, corresponding to a positive LCE 
in the secular evolution. In this paper, we refer to the stability for the system means that 
the orbiting planets can remain periodic or quasi-periodic motions with bounded trajecto- 
ries after the investigated time was done. We define the unstable orbits that either of the 
planets is ejected far away or moves too close to their parent star, and in our simulations, 
the integration was automatically ended when meeting the following criteria: (1) either of 
the eccentricities approaches unity, (2) either of the semi-major axes exceeds the factor of 2 
or reduces half of the starting values, (3) either of the planets collides with the star or these 
planets do each other when entering the scope of the mutual Hill sphere. 



3. Stable planetary geometry in 2:1 mean motion resonance 

To begin the study of the stable resonant geometry in HD 82943, let us firstly recall 
several known facts in our solar system. One may be aware that the pair of Jupiter and 
Saturn is in a near 5:2 commensurability. In general, the mean motion resonance takes place 
in the pairs of the moons of the major planets (such as the Galilean satellites of Jupiter: 
Io-Europa-Ganymede; see Lee & Peale 2002), the asteroidal belt and Kuiper Belt Objects 
(KBOs) (Duncan, Levison, & Budd 1995; Wan & Huang 2001). As more and more extrasolar 
planetary systems are being discovered, the resonant pairs are frequently found to occur in 
the multiple systems, e.g., HD 82943 (Gozdziewski & Maciejewski 2001), GJ 876 (Lee & 
Peale 2002), and possibly HD 160691 (Gozdziewski, Konacki, & Maciejewski 2003; Bois et 
al. 2003) in a 2:1 MMR, and 55 Cnc in a 3:1 MMR (Ji et al. 2003a), and so on. Several 
theoretical or numerical works have been done to enhance the comprehension of the evolution 
into the resonances by non- conservative tidal forces (Ferraz-Mello, Beauge, & Michtchenko 
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2003), the convergent migration for the planets due to the disk-planet interaction by a 
damping force (Snellgrove et al. 2001; Lee & Peale 2002; Kley 2003), the excitation of 
the orbital eccentricities by repeated crossings of two companions migrating on divergent 
orbits (Chiang et al. 2002; Chiang 2003) or the eccentricity growth by two resonant planets 
through inward migration (Murray, Paskowitz, & Holman 2002), and resonant inclination 
excitation mechanism for migrating planets in the gas disk (Thommes & Lissauer 2003). On 
the other hand, if the two planets are trapped into a 2:1 resonance via one of the potential 
mechanisms, it is still necessary to understand their actual situation of orbital motions after 
the capture, and to investigate whether they can secularly survive in the system in some 
appropriate geometry and examine how the stability depends on the orbital parameters. 

In the usual notation of celestial mechanics, the lowest order eccentricity-type (see 
Murray & Dermott 1999) resonant arguments 9 1: 9 2 for the 2:1 MMR are 

01 = Ai-2A 2 + tJ7i, (1) 

2 = Ai - 2A 2 + w 2 , (2) 

where Ai, X 2 are , respectively, the mean longitudes of the inner and outer planets, and wi, 
zu 2 denote their apsidal longitudes, respectively. Additionally, the relative apsidal longitudes 
of two companions #3 reads, 

9 3 = w\ — zu 2 = Aw. (3) 

As we mentioned previously, the apsidal alignment (or antialignment) is found to exist in 
most of the multiple systems (Kinoshita & Nakai 2000; Rivera & Lissauer 2000; Chiang, 
Tabachnik, & Tremaine 2001; Laughlin et al. 2002; Lee & Peale 2002, 2003; Gozdziewski et 
al. 2003; Ji et al. 2003a, b; Zhou & Sun 2003), to play a significant role in retaining the 
stability of the system (Gozdziewski 2003; Ji et al. 2003b). Moreover, this apsidal phase- 
locking of two orbits indicates that a pair of planets have common time-averaged rate of 
apsidal precession, and in the literature, it is also referred to apsidal resonance or corotation 
(Chiang et al. 2001; Malhotra 2002). For a 2:1 MMR, it is easily noticed that for the above 
three arguments, no more than two are linearly independent, indicating that either all three 
librate or only one librate in the occurrence of the libration of the critical argument. In more 
general case, Nelson & Papaloizou (2002) also showed that both situations can occur for 
a p : q eccentricity-type commensurability (Murray & Dermott 1999), where the associated 
resonant angles are defined by p , g ,A;(Ai, \ 2 , w±, w 2 ) = — q\ 2 + {k + {q — k)w 2 ){p^ q, k 

are positive integers and p < k < q) in the numerical surveys of the pair of the extrasolar 
planets. As a matter of fact, this was again verified in the case of the 3:1 MMR for 55 Cnc (Ji 
et al. 2003a; Zhou et al. 2004), where not only do all three resonant angles of A& — 3X c + 2w c , 
\b — 3A C + Wb + w c and A& — 3A C + 2wb (here the subscript b, c, respectively, denote the two 
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inner planets of 55 Cnc b and 55 Cnc c) librate, but the relative apsidal longitude of Azu 
librates about a constant angle simultaneously. 

In the simulations of HD 82943, we found that three types of stable orbits can survive in 
this system: (I) Only 9\ librates about 0°, while 9 2 and 9 3 circulate; (II) Case of alignment, 
three arguments all librate about 0°, denoting by 9\ 9 2 ~ 9% xs 0°; (III) Case of antialign- 
ment, both 9\ and #3 librate about 180°, while 9 2 librates about 0°, denoting by 6\ ~ 180°, 
9 2 ~ 0°, 9 3 ~ 180°. In the following sections, we will separately discuss above three stable 
configurations with both the numerical outcomes and analytical means and also compare 
them with remarkable evidences in other multiple systems involved in 2:1 resonance. 



3.1. Coplanar semi-analytical model 



First, let us suppose two planets with the masses of m\ and m 2 orbit the central star 
with the mass of m . Here we only take into account the point-mass interaction between 
star-planet and planet-planet, without the consideration of the effects of the oblateness 
and general relativity arising from the parent star. Thus, in Jacobi coordinates (cij,i = 
1,2; assume a± < a 2 ), the Hamiltonian for three-body system can be written (Brouwer & 
Clemence 1961): 

F = F + $ , (4) 



where 



and 



Gm mi Gm m 2 



2ai 



2a 2 



Gm\m 2 ^ /ll 

Gm m 2 

ri2 \ r 02 r 2 



(5) 



(6) 



Here r 02 (or r i2 ) is the distance between m (or mi) and m 2 , r 2 the distance of m 2 with 
respect to the center of mass for m Q and m\. Therefore, the disturbing potential $ can be 
further expanded to the following forms for the coplanar case, 



where 
Here 



Gm 1 m 2 ^ 

<£ = > Mai, a 2 , ei, e 2 ) cos0, 

a 2 ^ 



k\\\ + k 2 \ 2 + /c 3 G7i + k 4 w 2 , 



J2^ = o. 



(7) 

(8) 
(9) 
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In (7), S depends on the semi-major axes (in association with the Laplace coefficients) and 
the eccentricities of two planets and this expansion of the disturbing potential also requires 
the smaller values of (i=l,2) and the two orbits without intersection. In the case of GJ 
876, Lee & Peale (2002) expanded (7) and collected the terms up to e 3 , the truncation series 
may be valid for low values of e^, where e\ = 0.255 and e 2 = 0.035 (Laughlin & Chambers 
2001), without the occurrence of the crossing orbits. In a more recent work, Beauge & 
Michtchenko (2003) also developed an expansion method that can deal with higher values 
of the eccentricities up to a limit of e; ~ 0.5, but they suggested that it is still related 
to the situation of the proximity between the planets. However, as for HD 82943, the 
planetary eccentricities can be as high as 0.41-0.54 for Fit 1, or 0.18-0.38 for Fit 2, the 
development of (7) with respect to large eccentricities cannot be exactly applied, because 
the expanding series, consisting of a great number of the expansion items, either converge 
very slowly or become divergent in the collision circumstances. Therefore, herein we seek for 
a semi-analytical treatment (Kinoshita & Nakai 2002) for this question and then compare the 
analytical results with direct numerical integrations. We emphasize that this semi-numerical 
method corresponds to a first-order theory with respect to the small parameter e, where e 
represents the mass ratio of the Jupiter-like planet to its parent star and is normally of the 
order of ~ 10~ 3 in the planetary systems, and can apply to the planetary system with any 
eccentricities, even for the mutual crossing case. 

The Hamiltonian for the coplanar case is 



In order to keep the Hamiltonian form, we should use Jacobi coordinates to study this 
system. As the indirect part of the Hamiltonian F does not contribute to the secular part, 
we simply take the direct part. As the two orbits of HD 82943 are close to each other and 
may intersect (Ji et al. 2003c; Hadjidemetriou & Psychoyos 2003; Lee 2004), the ordinary 
analytical expansions of the main part invalidate in the practical usage, then we adopt the 
original form of the main part of the disturbing function and numerically evaluate it. Again, 
by eliminating short-periodic terms in the classical average process, the new Hamiltonian 
reads: 



The degree of freedom of the new Hamiltonian (11) is reduced from four to two. However, 
the Hamiltonian F* is not integrable. As we shall see in §3.2 that the semi-major axes of 
the two resonant planets slightly change when they are close to 2:1 MMR, so we can assume 
that a\ and a 2 are almost constant (see Figures 2 and 3) for the critical argument Q\ = 0° 
or 9\ = 180°. In final, the degree of freedom of the new Hamiltonian is reduced to one: 



F = F(a 1 , a 2 , ei, e 2 , w u w 2 , Xi, A 2 ). 



(10) 



F* = F*(a 1 ,a 2 ,e 1 ,e 2 ,6 1 ,e 3 ). 



(11) 



F* = F*(e 1 ,e 2 ,e 3 ). 



(12) 
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We again eliminate e 2 in equation (12) with conservation of the total angular momentum J 
(see also Appendix), then we have 

F* = F*( ei ,9 3 ,J). (13) 

Thus we can draw the level curves of the Hamiltonian and understand the global behavior 
of ei and 83. As the eccentricity e\ of the inner planet becomes large, we numerically 
averaged the original Hamiltonian (10) under the condition of the critical argument 8\ = 0° or 
Q\ = 180° and the angular momentum conservation relationship, and then obtain numerically 
the averaged Hamiltonian (13) with the parameter of the total angular momentum J. The 
technique of the numerical averaging is extensively described in the paper (Kinoshita & 
Nakai 1985), in which the case of Q\ 7^ is also discussed. To make a study of the evolution 
of the eccentricities and the relative periapse arguments, we can draw the contour map of 
Hamiltonian (13) by taking 83 as the horizontal axis and e\ or e 2 as the vertical axis (see 
Figure 4) with the parameter J, which is determined from the initial conditions. 



3.2. Results for Fit 1 

For the first 100 runs for Fit 1, each orbit integration lasts for t = 10 7 yr. In the 
simulations, we notice that all the stable cases are involved in the 2:1 MMR and easily 
understand that the stability of a system is sensitive to its initial planetary configuration. 
Furthermore, we noticed that 2%, 3%, and 2% of the cases belong to Types I, II, and III and 
93% of the systems destabilize, which means that either of the two planets leave in situ to the 
infinity with a rapid increase of the semi-major axes or the eccentricities grow much larger 
on the timescale ~ 10 4 — 10 5 yr or even shorter, owing to the mutual interaction between 
star-planet or planet-planet through frequent close approaches. In Table 1, are listed two 
sets of orbital parameters-Fit la is related to the antialigned configuration for Type III and 
Fit lb is involved in the case for aligned orbits for Type II. Here the adopted parameters 
for two planets are the masses of 1.63-0.88 M jup , the semi-major axes of 1.16-0.73 AU, the 
eccentricities of 0.41-0.54. Let us first examine the stable orbits of Type I, where in the 
numerical simulations, we find that only one argument of 8\ librates about 0° and the others 
82 and 83 take up the full circulation from 0° to 360°. The semi-major axes a± and a-i perform 
tiny oscillations about 0.73 and 1.15 AU for the whole integration time. Type I was initially 
introduced by Gozdziewski & Maciejewski (2001) with the previous fits. However, in this 
paper, Types II and III, which can grasp the couple of companions escaping from deviating 
their tracks in two directions of remaining both the semi-major axes and the eccentricities, 
are new findings of stable orbits for this resonant system. Hence, in the following sections, 
we primarily pay attention to the stable configurations linked to both a 2:1 resonance and 
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the apsidal corotation. 

3.2.1. Antialigned case 

Figure 2 shows the typical orbital evolution for Fit la (antialigned orbits). In the figure, 
we can see that the semi-major axes a x and a 2 slightly vibrate about 0.73 and 1.15 AU for 
t = 10 7 yr 5 (in the figure, to see more clearly, we simply display a snapshot for t = 2000 
yr), meanwhile 8\ librates about 180° with a moderate amplitude of ~ 30°, 8 2 about 0° with 
a smaller amplitude of ~ 10° and 63 (by thick line) about 180° with an amplitude of ~ 30°. 
It is noteworthy to point out that the 2:1 resonance variable #1 librates in accord with the 
relative apsidal longitudes 83, occupying the same period of ~ 300 yr. This phenomenon can 
be easily understood because the above three angles are not independent in the course of the 
orbital motion. Thus we can call 9\ and 83 as two fundamental variables to characterize the 
stable planetary geometry in 2:1 MMR. In addition, the modulations of the semi-major axes 
a; and 9i (i = 1,2) reveal the fact that the two planets are indeed trapped in a 2:1 MMR 
for the secular orbital evolution. In Figure 2, we further notice that the eccentricities of e\ 
and e 2 separately range in (0.36, 0.56) and (0.39, 0.49), bearing the equal libration periods 
of 81 and 83. The orbital antialignment of axes in HD 82943 is reminiscent of that in the 2: 
1 orbital resonances for Io-Europa system (Lee & Peale 2002), in which 8\ (involved in the 
Io's longitude of periapse) librates about 0°, and 8 2 (involved in the Europa's longitude of 
periapse , #3 librate about 180°, respectively. The antialigned circumstance for Io-Europa 
is equivalent to Type III in the sense of dynamically stable configuration, however, these 
jovian moons are occupying almost near-circular orbits resulting from the primordial or the 
tidal origin (Peale 2003) and the resonance configuration may come into reality when the 
eccentricities are small. By contrast, the HD 82943 system differs the Io-Europa pair in 
that both of the planets possess high eccentricities at present day. Hence, the origin of the 
2:1 resonance for HD 82943 (Fit la) orbiting two massive planets with elliptic trajectories 
should be explained by an innovative mechanism rather than a capture induced by co-orbital 
differential migration (S. Peale 2003, private communication; see also §5) and perhaps stem 
from gravitational scattering in a crowded system that augment their eccentricities, i.e., for 
a system with 10 planets (Adams & Laughlin 2003), the scattering processes can yield the 
full extent of possible eccentricity where < e < 1. 

Recently, Bois et al. (2003) found that the stable configuration of the HD 160691 



5 We also extend the integration to 10 8 yr using symplectic integrator for this case, and no sign shows 
that the system will become chaotic for longer evolution. 
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system is possibly involved in a 2:1 mean motion resonance combined with an anti-aligned 
in two apsidal lines, with higher eccentricity (e 2 > 0.52) of the outer planet. Obviously, their 
research supported the stable geometry of Type III of HD 82943, which also implies that this 
resonant topology may be suitable especially for the planets bearing larger eccentricities. 

Moreover, this antialigned resonant configuration for Fit la means that the conjunctions 
can take place when the inner planet moves near apoapse and the outer planet is close to the 
periapse, and vice versa, which indicates that the two planets are far from each other during 
their orbital evolution and protected from frequent close encounters. Therefore, we can safely 
conclude that the stability of the HD 82943 system, which is related to the above coplanar 
configuration for Fit la, can be simultaneously sustained by two dynamical mechanisms: the 
2:1 MMR and the apsidal antialignment. 

3.2.2. Aligned case 

Figure 3 exhibits the orbital evolution for Fit lb (aligned orbits). Consequently, let 
us note that the semi- major axes a\ and a 2 do not change dramatically but undergo small 
oscillations about 0.73 AU and 1.15 AU for t = 10 7 yr (a snapshot for t = 3000 yr), further 
the eccentricities e\ and e 2 librate in the extent of (0.5, 0.8) and (0, 0.45), respectively, in 
the behavior of the converse cycles, owing to the conservation of the total angular momenta 
for the system. Here Q\ librates about 0° with a small amplitude of ~ 20°, but 6> 2 and 6* 3 (by 
thick line) individually librate about 0° with a large amplitude of ~ 70°. In this resonant 
aligned geometry, we again observe that 9±, 9 2 and 9 3 share the common librating period of 
~ 700 yr, in the same time, coupled with the librating period of the eccentricities of e\ and 
e 2 . In addition , it is worthwhile to point out that at the time that 9±, 9 2 and 9 3 all approach 
zero (where 9\ and 9% start from the negative values to zero), the eccentricity of e\ reaches 
the maximum ~ 0.8 and e 2 is close to the minimum ~ 0, thus at the aligned conjunction 
point, the separation between two planets d = a 2 (l — e 2 ) — ai(l — e x ) (where a± and a 2 are 
almost preserved on the condition of 2:1 resonance locking) can reach the maximal value 
and the maximum separation can prevent them from strong mutual perturbation when they 
rendezvous at the synodic location. Similar analysis can be applied for the case of Fit la, 
where both 9\ and 9% run through 180° (see Figure 2). Again, let us make a brief comparison 
with the circumstance for the 2:1 resonant system GJ 876, which is also found to be in aligned 
configuration for the two orbits (Marcy et al. 2001). Lee & Peale (2002) further found that 
the amplitude of the librations of three arguments are not far from 0° by using Laughlin- 
Chambers solutions and they even indicated that all the librations are remained even for 
the amplitudes of 9\ amounting up to 45° for some cases. As for Fit lb, let us bear in mind 
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the fact that the planetary eccentricities are much larger than those for GJ 876, although 
they may have similar planetary orbital geometry Sequently, the conjunctions always take 
place when two planets move close to their the periapse longitudes and the alignment can be 
maintained over long-term orbital evolution. This, de facto, reveals that Type II coplanar 
configuration may be possible for real observed systems of GJ 876 and HD 82943. 

3.2.3. Comparison with the semi- analytical results 

Here we adopted the semi-analytical model introduced in §3.1 to draw the Hamiltonian 
contour and then compare these results with those given by direct integrations over secular 
orbital evolution. In the averaged Hamiltonian after removing the short-periodic terms, the 
difference arising from the choice between the Jacobi coordinates and astrocentric coordinates 
is the order of the perturbation (mi /mo, m 2 /m ), which could be ignorable. Therefore, we 
adopt the astrocentric elements given in Table 1 to be Jacobi elements in the averaged model. 
As a comparison, the initial orbital elements for numerical integration are also converted into 
the Jacobi coordinates. In Figure 4, we indicate that the contour diagram of Hamiltonian 
(13) are plotted against various levels by changing e±, e 2 and Aw = zui — w 2 , left panel 
represents the results for Fit la, where the contour levels are exhibited by thin line and 
the numerical solutions shown by thick line (the innermost curve about marked plus), while 
right panel for Fit lb. The — Aw figures show a good agreement between the numerical 
outcomes and the semi-analytical solutions for two fits. It is worthwhile to mention that the 
closest level of the Hamiltonian contour to the numerical solution set the stable boundary 
for the apsidal resonant geometry, inside which the regular orbital motion for two orbits 
can be found with qualified eccentricities and the relative apsidal longitudes, thus we may 
conclude that this semi-analytical means is effective to help predict stable orbital solutions 
(e.g., Beauge et al. 2003). On the other hand, the diagrams of the Hamiltonian contour 
further show that the stability zones for the eccentricities ei and e 2 are connected with the 
librating amplitude of Aw: the smaller stability zones of the eccentricities, the narrower 
modulations of the relative apsidal longitudes, and vice visa. 

In addition, we calculated the equilibrium points for each fit: (180°, 0.45) and (180°, 
0.46), respectively, for Fit la; (0°, 0.68) and (0°, 0.29), respectively, for Fit lb. As for 
Fit lb, we observe that there should exist stable orbits about the equilibrium center for 
(Aw, e 2 ) = (0°, 0.29), together with high eccentricity ei above 0.60, which satisfy the neces- 
sary conditions for two planets evolving into a 2:1 resonance (Kley, Peitz, & Bryden 2004; 
see also §5). Again, we recall that the period of the 2:1 resonant variable 0i is equal to that 
of the relative apsidal longitudes 9 3 (see Figs. 2 and 3), which can be clearly reflected in 
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Fig. 4. Moreover, in a view of the numerical results, we can see that the eccentricities of 
both planets are well restricted with the small-amplitude libration of 3 about 180° for Fit 
la (antialigned) and large-amplitude libration of 9 3 about 0° for Fit lb (aligned), which can 
also be consistent with the results in Figs. 2 and 3. 

As aforementioned, different selections for the initial orbital parameters will lead to 
stable or unstable configurations for the HD 82943 system. By analyzing the starting orbital 
data, we discover that those stable planetary orbits involved in Type II (or III) require the 
initial 9i and 8% should be constructed to satisfy the conditions of being less than tens of 
degrees from 0° (or 180°), implying that the linear combinations of the mean anomalies, 
ascending nodes and periapse arguments can be close to the 2:1 resonant geometry. 

3.3. Results for Fit 2 

As for the second 100 runs near the vicinity of Fit 2, we found there exist two kinds of 
stable geometry related to a 2:1 MMR for this new fit (Mayor et al. 2004), in which 16% 
of the stable systems can last for 10 7 yr (where 14% of the stable cases is in aligned orbits 
of Type II and 2% with simply 6i libration for Type I), and 84% are found to self-destruct 
at t < 10 6 yr. The percentage of the survival orbits indicates that there is a likelihood for 
the two planets of this system to be aligned. In Table 2, are listed the orbital parameters 
in association with the alignment of two orbits. Let us recall the adopted planetary masses 
of 1.84-1.85 Mj up , the semi-major axes of 1.18-0.75 AU and the eccentricities of 0.18-0.38. 
In the following subsections, the additional computations are carried out most frequently 
according to the initial orbits from Fit 2, and the typical secular evolution timescale ranges 
from 10 6 yr to 10 7 yr. 

Figure 5 illustrates the typical aligned orbital evolution for Fit 2. In comparison with 
Figure 3, we notice that the semi-major axes a\ and a 2 are almost unchanged but modulate 
about 0.75 and 1.18 AU with small amplitude for t = 10 7 yr, at this time, the amplitude 
of the oscillations for e\ and e 2 are not so large and the eccentricities are just wandering in 
the span (0.34, 0.44) and (0, 0.25), respectively. Here Q\ librates about 0° with a moderate 
amplitude of ~ 45° (in coincidence with one of the aligned case for GJ 876 found by Lee 
& Peale 2002; see also §3.2.2), but 9 2 and 9 3 (by thick line) individually librate about 0° 
with large amplitude of ~ 80°. Again, we notice that the eccentricities have the libration 
periods of ~ 600 yr, which are coupled with those of 6\ and #3. It should mention that 
although the various starting orbital parameters between Fit 2 and Fit lb can result in 
variational librating amplitudes in the eccentricities and resonant arguments, however, from 
a viewpoint of topology, there is not so much difference between two aligned orbital solutions. 
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Furthermore, in Figure 6, we still plot t{ — Azu Hamiltonian contour (by thin lines) and the 
numerical results are presented by dotted points. The equilibrium solutions (by plus) are: 
(0°, 0.40) and (0°, 0.13), respectively. The semi-analytical contour chart exhibits a good 
agreement with the numerical investigations of maintaining the eccentricities with larger 
#3-libration amplitudes. On the other hand, the figure implies that the ^3-libration will be 
broken up if the libration amplitudes exceed or approach the critical value of 90° near the 
separatrix, then two successive scenarios may occur, either aligned orbits evolve into Type 
I geometry that simply 61 librates, or the system destabilize on the shorter timescale. In 
the sense of the present solutions, the aligned configuration for HD 82943 is not robust than 
the antialigned case, where a wider relative periapse longitudes can be initially adopted (see 
Figure 4). Bois (2003, private communication) also confirmed that the stable strip of the 
alignment for the initial W\ — w 2 is much narrower than that of the antialignment in their 
numerical calculations for HD 82943. However, if the aligned orbits are close to librating 
center, the robustness can be strengthened with a different story. 

In order to distinguish the dynamical behaviors in orbital parameter spaces and explore 
how the stability zones dependence on these parameters, we extensively perform abundant 
integrations for each series (where in each series, we only allow one parameter vary but others 
keep unchanged) to numerically search for the stable solutions. In the following sections, we 
will discuss these results respectively. 

3.3.1. Dependence on the semi-major axes 

Figure 7a show the numerical scanning of the HD 82943 system in the semi-major axes 
space of [ai,ai]. All other orbital parameters are taken from Table 2. In this group of 
the simulations, we choose a x G [0.60 AU, 0.80 AU] and a 2 G [1.00 AU, 1.30 AU] with the 
resolution of 20 x 20 grids, each direction spaced 0.01 AU. In Fig. 7a, the cross indicates the 
unstable orbits and the filled circles for stable orbits that the planets survive at least 10 6 yr. 
Here, we notice that there exists a stable strip with the width from ~ 0.01 AU to ~ 0.02 
AU, indicating that these stable solutions are quite close to the 2:1 resonance. Meantime, 
the intersection between the vertical and horizontal lines represents for the present solution 
from Fit 2, which is well located in the stable map. 

In the case of 2:1 resonance, we have a 2 ja\ = (P 2 /Pi) 2 ^, where P 2 jP\ ^ 2, if both of 
orbits contract or expand k times of the initial configuration (e.g., Fit 2) and the planetary 
masses becomes y/k of the initial masses (where other orbital parameters remained), the 
semi-amplitude of the radial velocity will not change, that's the reason why we can see 
a linear stable ribbon in Fig. 7a. In this sense, the semi-major axes cannot be uniquely 
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determined unless their orbital periods are previously obtained from the observations. 

3.3.2. Dependence on the eccentricities 

Mayor et al. (2004) pointed out that one should be careful of employing the newly 
updated orbital fits (Fit 2) due to the derived best-fit solutions of lack of considering the 
planet-planet interaction. Hence, the exploration for the orbital parameter space is extremely 
useful to help understand the stability regions relying on the eccentricities. Similar to the 
numerical scanning in the semi-major axes space, here only both of the eccentricities of e± 
and e 2 are allowed to be free parameters. In this series of the investigations, we let e± G 
(0.0, 0.60 ) and e 2 G (0.0, 0.60) and the resolution in [ei, e 2 ] space is 20 x 20 points, spaced 
0.03 in each grid (see Figure 7b), the integration time is 10 6 yr for each orbit. In the figure, 
the symbols are defined as in Figure 7a, where the filled circles mean those orbits are steady 
over secular evolution. 

According to the eccentricity of the outer planet e 2 , we summarize the results as follows: 

(1) strongly stable cases, where < e 2 < 0.24 (40%), 

(2) partially stable cases, where 0.24 < e 2 < 0.39 (16% stable, 9% unstable), 

(3) unstable cases, where 0.39 < e 2 < 0.60 (35%). 

As for case (1), we should mention two points: firstly, the eccentricities of Fit 2 (shown 
by the crossing of two lines) just reside in the stable map of [e±, e 2 ]; secondly, in this figure, 
we notice that there exist stable solutions when e\ = 0.4 and e 2 ~ 0, which is consistent with 
Mayor et al. (2004), who also provided the other solution involved in an aligned configuration 
with a near-circular orbit for outer planet but found the residuals do not dramatically change. 
We can also note that a higher e\ can be possible in this system. Both of cases (1) and (2) 
contain more that 50% of the stable orbits in the numerical scanning, exhibiting that only a 
moderate e 2 will favor the stability of the HD 82943. 

The case (3) shows that when e 2 becomes larger than 0.39, no regular orbits will appear 
even though e^ can be adopted to be small values, which suggests that in this situation the 
stability for HD 82943 is more sensitive to the orbital motion of the outer planet, which can 
greatly determine whether the system is in chaotic or regular status. 
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3.3.3. Dependence on the mean anomalies and arguments of periapse 

Figures 7c and 7d exhibit our numerical exploration in the parameter spaces of [Mi, M 2 \ 
and [ui, u 2 ] for the time span of 10 6 yr, respectively, with the resolution of 18 x 18 grids. The 
stable map in the [Mi, M 2 ] space again show the condition of the 2:1 MMR , where there is 
a linear relationship M x = 2M 2 + const. Therefore, on the basis of Fig. 7c, we can not only 
verify our original mean anomalies of Fit 2 (by intersection) but find the other arguments 
related to stable orbital motions for two planets . The stable diagram in the [cji,cj 2 ] space 
reveals the fact that both of the apsidal longitudes of w\ and w 2 (where nodal longitudes 
are constant here) precess at the equal rate about the star. The stable zones in Fig. 7d 
are quite narrower, as we stated before, this results from the aligned configuration with a 
larger libration of 63 ~ 80° (ref. Fig. 5). Both figures clearly illustrate that those stable 
configurations are linked to the 2:1 resonance and apsidal corotation. 

3.3.4- Dependence on the planetary mass ratio 

In this section, we expect to understand how the stability depends on the planetary 
mass ratio for the coplanar planetary configurations. Here in our numerical study, we still 
adopted the orbital elements from Table 2, except the masses of two planets mi and mi. 
Firstly, we varied sini in an incremental step of 0.05 ranging from 0.35 to 1.00, and the 
rescaled masses of m k (k — 1,2) can conveniently obtained by the product of m^suii for a 
specific sini, where the mass ratio m 1 /m 2 is fixed to be 1.85/1.84 ~ I. Then 14 coplanar 
integrations were carried out for 10 7 yr in this runs. The numerical results show that the 
stable orbits (of Type II) for HD 82943 requires sini > 0.50. The stability set limitations on 
the planetary masses where mi e [1.85Mj up , 3.70M/ up ] and m 2 e [1.84M/ up , 3.68M/ up ] for 
different values of sini. In the case of GJ 876, similarly, Marcy et al. (2001) reported that 
several stable systems, catalogued in Type II, had sini = 0.50 for the integration timescale 
up to 500 Myr. Additionally, Laughlin & Chambers (2001) also suggested that sini possibly 
ranges from 0.5 to 0.8 for the coplanar fits for GJ 876. From above statements, we notice 
that sini > 0.50 is necessary for both systems occupying a 2:1 resonance and the extent of 
the planetary masses can further be estimated. 

In the other computations, we changed the mass ratio for m\jm 2 and fixed the mass 
m 2 = 1.84 Mj up , let m\/m 2 G [1/4,4], then performed over thirty additional simulations to 
explore the stability of the system depending on the variational mass ratios. As a result, we 
found that m 1 /m 2 < 2 is necessitated for the stable orbits, indicating an upper limitation for 
the mass of the inner planet when we remain that of the outer companion. In the following 
step, we replaced mi with much lower masses ranging from several M e to 20 M e for new 
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investigations, and found that the stable geometry of Type II can still be remained and a 
surprise is that the eccentricity of the inner planet can amount up to ~ 0.90 with quasi- 
periodic oscillations. These numerical explorations suggest that the apsidal configuration can 
maintain with much smaller masses for inner planet, which is consistent with the results given 
by Beauge et al. (2003), who analyzed the corotational solutions by gradually decreasing 
the mass of the inner planet and found that although the solutions simply depend weakly 
on the individual masses of the planets and the maximum value for the planetary mass of 
mi is indeed in existence with given eccentricities for the stable apsidal topology. 



3.3.5. Dependence on the relative inclination 

In our solar system, most of the major planets have small inclinations with regard to 
the essential reference plane. However, there is a special case for the Neptune-Pluto system 
(M Neptune pa ^M Saturn ), where Pluto orbit an eccentric trajectory with a large inclination 
of ~ 17° but Neptune travels on a near-circular orbit with the inclination less than 2° 
(Murray & Dermott 1999). Cohen & Hubbard (1965) discovered that these two planets are 
involved in a 3:2 resonance with the argument 3Ap — 2Aat — (^p,n denote, respectively, the 
mean longitudes of Pluto and Neptune, wp is the longitude of perihelion of Pluto). These 
studies reveal that there may exist a relative inclination configuration between two orbits in 
association with the mean motion resonance. In this sense, we anticipate to understand the 
situation for the present exosystem. 

We again performed extra tests of mutually inclined orbits for two resonant planets for 
the integration of 10 7 yr, to examine the stability of the HD 82943 system. The orbital data 
were taken from Table 2 and for simplicity, we just adopt their minimum masses and remain 
them unchanged for our additional integrations. For the initials, here we simply altered the 
inclination of the outer planet li and increased I2 in an increment of 5° from 5.5° to 85.5°, but 
always assumed the inclination of the inner planet l\ to be 0.5°. The numerical integrations 
reveal the fact that although the inclination was greatly changed, the two planets are still in 
2:1 MMR and apsidal alignment (Type II) with the inclination of I 2 < 25.5°. On the contrary, 
for I2 > 25.5°, we found that the eccentricity of the inner planet can pump up to unity due 
to the great perturbation exerted by the outer planet and leave the orbits eventually. The 
relative inclination i r can be determined by cosi r = cos/icos/2 + sin fx sin/ 2 cos(f2i — f^), 
where Q\ and Sl 2 are, respectively, the longitudes of ascending nodes for two planets. To 
sum up, in the case of the non-coplanar configuration for HD 82943, stability requires that 
the relative inclination be ~ 25° or less. In the case of GJ 876, Rivera & Lissauer (2001) 
carried out dynamical fits for the possible best solutions for the mutual inclination geometry 
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and found that the stability for GJ 876 requires the relative inclination to be 12° or less. 
Therefore, it should be mentioned that the mutually inclined orbits dynamically survive over 
long-term timescale in the planetary geometry locking into a 2:1 MMR, but it should place 
constraints on the initial conditions for two planets. Thommes & Lissauer (2003) showed that 
two planets migrating in resonance can perform significant excitation of their inclinations and 
the relative inclination between two planets can grow to ~ 30° in their three-dimensional 
simulations. As a final note for this part, the above-mentioned numerical investigations 
for various mutual inclined exosystems show that the low relative inclination do play an 
important role in the stability for such systems, and this mutual inclined configuration may 
arise from the inclination-type resonance (Thommes & Lissauer 2003) associated with the 
longitudes of the ascending nodes for inner and outer planets if the inner planet owns a high 
eccentricity (e.g., for Fit 2, e\ = 0.38), while the higher relative inclination geometry can 
destabilize the system (e.g., Lissauer & Rivera 2001) on shorter timescales. 



3.4. Comparison with other works 

As we mentioned previously, several studies have been devoted to explore the orbital 
motions of the exosystems in 2:1 MMR. In the case GJ 876, Lee & Peale (2002) discovered 
that there are stable configurations with 9i, 9 2 and 83 all librating about 0° (related to Type 
II) for 0.15 < e 1 < 0.86, and Kinoshita & Nakai (2001) and Ji et al. (2002) also studied 
this system and confirmed the 2:1 resonance acting as an effective mechanism to hold the 
system, moreover, an extended stable geometry of Type I was shown by Ji et al. (2003b) 
for GJ 876. In other works, Hadjidemetriou & Psychoyos (2003) numerically studied the 
families of periodic orbits for HD 82943 in the rotating frame, and indicated that two kinds 
of families orbits can survive for the 2:1 resonant planets in their computations. However, 
it is the first time that we show an exhaustive exploration in the orbital parameter space 
for the HD 82943 system based on the updated solutions, which reveals the new dynamical 
features for this system. 

In summary, the results presented in this article on the possible stable planetary geom- 
etry involved in a 2:1 resonance are not only consistent with the previous studies, but also 
outline an extensive dynamics of these systems. Moreover, we point out that the aforemen- 
tioned three groups of resonant configurations imply a general principle that governs two 
planets in a dynamically stable state over secular orbital evolution in the planetary systems, 
which should be examined by more future discoveries of the 2:1 resonant pairs. Still, we 
will concentrate on whether a postulated terrestrial can survive in such geometry (see §4), 
and how a certain system is evolved into an apsidal alignment or antialignment resonant 
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configuration in §5. 

4. Habitable zones 

The Habitable zones (HZ) are generally convinced to be suitable places for terrestrial 
planets that can provide the liquid- water, subtle temperature and atmosphere environment, 
and other proper conditions (Kasting et al. 1993), supporting the development and biological 
evolution of life on their surfaces. The habitable region is considered to be centered at a p 
= 1 AU (L/Lq) 1 / 2 (Gould, Ford, & Fischer 2003) or a p > 1.03 AU (L/L Q ) 1 ' 2 (Kuchner 
2003), such that a water planet can exist in equilibrium with stellar radiation, where a p is 
the radius of the planet's orbit and L is the luminosity of the host star. In a viewpoint of 
the dynamics, first of all, there should be stable orbits at ~ 1 AU for Earth-like planets 
moving about main sequence stars over very long timescale and we can also refer this to the 
dynamical habitability. In a recent work, Menou & Tabachnik (2003) exhaustively studied 
the dynamical habitability of 85 known extrasolar planetary systems through numerical 
simulations of their orbital dynamics in the presence of potentially habitable terrestrial 
planets. They found that more than half of the known exosystems cannot harbor habitable 
terrestrial planets and about 25% of the systems mostly with close-in giant planets are 
dynamically habitable, similar to our solar system. And in this section, we only draw our 
attention to the potential habitable zones for the systems (e.g. HD 82943 and GJ 876) with 
a stable geometry in a 2:1 resonance. 

Here we performed numerical surveys in examining the Habitable zones both for HD 
82943 and GJ 876. And at first, we considered stable configurations for HD 82943 (Fit 2) 
and GJ 876 (Keck and Lick Fit with sini = 0.78; see Table 2 in Laughlin & Chambers 
2001), which can be viewed as genuine systems closest to the actual observations. Then, 
we generated 100 seed planets that all bear the same masses of Earth (where M e ~ 3.14 x 
lO~ 3 M jup ) in the simulations for each system. In Figure 8, we can see the distribution of the 
initial orbital elements for these postulated planets-they initially move on the less inclined 
belt with the relative inclination with regard to the plane of the resonant pair no more than 
5°, the semi-major axis a G [0.96 AU, 1.05 AU], the eccentricity e e [0,0.1], and the leftover 
arguments are set at random. Here we simply accounted for the near-circular orbits rather 
than much more eccentric orbits, because the life should develop and evolve in the biological 
environment without larger variations of temperature, implying the Earth-like planet should 
not move neither too close nor too far from the parent star. 

For each star- two-planet-" Earth" four-body system, we carried out the integration for 
the timescale of 10 7 yr, to examine the dynamical habitability for the HD 82943 system. 
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As we shall see, we find that none of the orbits can retained stable in the region at ~ 1 
AU for the integration time, the seed planet travels on hyperbolic or parabolic orbits at 
the characteristic timescale r < 2 x 10 5 yr (where 81% of them are ejected at r < 5 x 10 4 
yr and 95% at r < 1.0 x 10 5 yr, see Figure 9a), owing to the scattering (Lin & Ida 1997) 
between two close resonant planets for the excitation in the eccentricity. A typical orbital 
evolution for the ejected orbits is shown in Figure 9b, the semi-major axe a grows from 
~ 1 AU to over 100 AU, the eccentricity e undergoes a rapid increase from ~ 0.1 to 1 and 
the inclination i is also excited to a high value of ~ 40° (in some cases, the inclination 
can exceed 90°), finally, at 4,800 yr, the assumed Earth-like planet is scattered away from 
the two-planet system, while the orbital motions of the resonant planets remain as usual. 
The gravitational scattering on the terrestrial planet caused an unstable motion either as an 
ejection or a catastrophic collision with the star or planets due to the dynamical instabilities 
(Ford, Havlickova, & Rasio 2001). Our numerical investigations suggest that there will be 
least possibility for stable Earth-like orbits surviving in the range of [0.73 AU, 1.18 AU], even 
for wider habitable regions [0.75 AU, 1.40 AU] given by Kasting et al. (1993), which is in 
good agreement with the simulations by Menou & Tabachnik (2003). Analytically, Gladman 
(1993) attained the minimum separation A between two planets moving on initially circular 
orbits require for stability, 

A = 2V3R H , (14) 

where the mutual Hill radius Rh is 

f mi + m 2 \ 1/3 fa x + a 2 \ 

Rh -{-^t) {-^) 

For Fit 2, we achieve A ~ 0.34 AU and Rh — 0.10 AU, which these values may become a 
bit different as the two planets for HD 82943 orbit on eccentric paths. Nevertheless, our goal 
here is to present a qualitative analysis, and we can easily notice that the initial orbits for the 
assumed planets are quite close to Hill sphere regime, resulting in their unstable motions. But 
it cannot rule out other stable regions for additional companions. For example, Sandquist 
et al. (2002) explored a three-planet case for HD 82943, which a supposed jovian planet 
with the mass of 0.5 M jup at 0.02 AU was added to the system, and they found that this 
additional planet does not affect the orbital motion of the two outer planets and also not 
disturb the stability of this system. On the other hand, for the exterior domain extending 
to 30 AU ~ 50 AU or even farther out, does there exist the dust rings or circumstellar disk 
that are analogous to the Kuiper Belt being a remnant of the solar nebula? The recently 
launched SIRTF mission will provide fruitful information on this point in its ongoing tasks. 

As a comparison, we accomplished the other same runs for GJ 876. And we conducted 
the four-body integrations for the time span of integration of 10 6 yr (~ 10 7 orbital periods 
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for the innermost planet) and found all the systems are stable for the integration timescale. 
Figure 10 exhibits the typical orbital evolution for the Earth-like planet, both the semi-major 
axis a and the eccentricity e execute small fluctuations about 1 AU and 0.06, respectively, 
and the inclination % also remains less than 2 degrees in the same time span. And there is 
no sign to indicate that such regular orbits at ~ 1 AU with low eccentricities will become 
chaotic for much longer time even for the age of the star. Due to the faint luminosity for 
M4 star GJ 876, the habitable zones may reside in the range of [0.10 AU, 0.20 AU] (Kasting 
et al. 1993; Menou & Tabachnik 2003) sufficiently nearer to the orbits of GJ 876 b (~ 0.13 
AU) and GJ 876 c (~ 0.21 AU) (Marcy et al. 2001; Laughlin & Chambers 2001; Rivera & 
Lissauer 2001), hence, the global chaos, corresponding to the circumstance of a collisional 
disruption or ejection, can always take place according to the Hill stability criterion (Gladman 
1993; Chambers, Wetherill & Boss 1996) or direct numerical examinations, and exclude the 
possible existence of Earth-like planets for such habitable regions. Whereas Kuchner (2003) 
pointed out that a liquid planet with a p > 1.03 AU (L/L ) 1//2 (where for GJ 876, a p > 0.12 
AU) may exist under the protection of a thick steam atmosphere, it seems plausible for 
the Earth-like planets revolve around GJ 876 at the orbit about 1 AU, in the concept of 
the dynamical habitability. However, our predictions should be confirmed through careful 
investigation in future research if the Earth-like planets do survive secularly in such habitable 
zone. Firstly, the direct evidence of the presence for such planets strongly depend upon the 
future high resolution astrometric measurements , e.g., Space Interferometer Mission (SIM) 
or Terrestrial Planet Finder (TPF) that can provide much better precision than the present 
long-term precision of ~ 2-3 m s" 1 (Butler et al. 2003) for ground-based Doppler surveys. 
Secondly, the simulations of the runaway accretion of the planetary embryos (Laughlin et 
al. 2002) suggest that it is impossible for the Earth-mass planets to be formed in Habitable 
zones even for a well-separated system with near-circular orbits (e.g., 47 Uma), because 
the formation of the terrestrial planets are terrifically constrained by giant planet migration 
(Armitage 2003), therefore the depletion of small planetesimals in a disk-clearing phase 
may possibly take place when the massive planets sweep over them in their orbital decay 
process unless this orbital migration can halt at a proper position allowing for the continuous 
growth of the planetary embryos or the already-created Earth-like planets can also follow the 
migration with the giant planets in step and cease about the Habitable regions. Nevertheless, 
fortunately, the terrestrial planets are conceivable to survive in a modest fraction of systems 
where a single generation of massive planets formed (Armitage 2003) without significant 
migration. 
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5. The Possible Origin for the Orbital Geometry 

The short-period giant planets that are close to the parent star with high temperature 
(e.g., 51 Peg, so-called hot Jupiter with a < 0.1 AU), are now thought to be formed through 
the shrinking orbital migration rather than be formed in situ due to the tidal interaction of 
a protoplanet with a surrounding gaseous disk (Ward 1997) and this orbital migration can 
originate from the exchange of angular momentum between the natal disk and protoplanet 
(Papaloizou 2003 and references therein). Now, it is widely believed that the migration will 
ubiquitously come about when two-planet systems suffer tidal interaction with an exterior 
disk, the outer planet will migrate inward until it gets trapped in a low order mean motion 
resonance with the inner companion (Bryden et al. 2000; Kley 2000) and the eccentricities 
will stop growing and be balanced as steady equilibrium values (Nelson & Papaloizou 2002) 
after being captured into the resonances. 

For the aligned orbits of Type II, in the earlier works for GJ 876, Snellgrove et al. (2001) 
performed full hydrodynamical two-dimensional simulations to study the inward migration 
for GJ 876 planets interacting with each other and a protoplanet ary disk, and the outer planet 
can be trapped into a 2:1 commensurability with the inner planet through slow migration, 
where the orbits were observed to be aligned with the resonant angles librating about zero 
(see also an overview in Papaloizou 2003) and the system can survive at least 2 x 10 7 
orbits by the removal of the external disk. Subsequently, Lee & Peale (2002) showed that 
the forced inward migration of the outer companion of the GJ 876 system leads to certain 
capture into the observed resonances if the initial eccentricities are smaller enough where 
ei < 0.06 and e 2 < 0.03 with proper migration rates. They also revealed that after resonance 
capture, the eccentricities increase rapidly due to resonant interactions between the planets 
and the forced migration without eccentricity damping, while with parameterized eccentricity 
damping where e^/e^ = —K\ai/di\, the eccentricities can reach equilibrium values that remain 
unchanged for sufficiently long migration in the resonances. Recently, in the study of the 
evolution in resonant planetary systems (e.g., HD 82943), Kley et al. (2004) found that the 
two planets can enter into a 2:1 resonance in varying timescale provided that e 2 < 0.25, 
with both full hydrodynamical simulations and damped N-body calculations. They further 
suggested that there are several requirements for the 2:1 resonance capture-the more massive 
mass for the outer companion, the higher eccentricity of the inner planet and the apsidal 
alignment of 83 = 0°. In the case of Type II (Fit 2), where the eccentricity e 2 = 0.18 and 
a moderate t\ = 0.38 and with two nearly equal planetary masses, or in the case of Fit 
lb, where there can also be stable geometry about e 2 ~ 0.25 and a higher e± with a larger 
outer planet, in both cases the two planets undergo # 3 -libration about 0°, thus we may 
comment that if the HD 82943 system is most likely to be aligned (Lee M.H. 2003, private 
communication), then the origin for such orbital geometry can be possibly attained by the 
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predictions of the hydro dynamic models (Kley et al. 2004), we will address this problem in 
our forthcoming paper. 

For the antialigned orbits of Type III, Thommes & Lissauer (2003) argued that if the 
planets reached this configuration by migration in resonance, a significant mutual inclination 
between their orbits would seem to be required. As in this geometry, the two orbits for HD 
82943 can experience the intersection on the co-orbital region (Ji et al. 2003c; Lee 2004), 
therefore a moderate inclination between two orbits can be achieved to act as a protection 
mechanism against close approaches (Thommes & Lissauer 2003). In §3.3.2, we also show 
that the highly inclined resonant configuration for this system can also be in presence, which 
the system is stabilized by both the 2:1 resonance and the apsidal phase-locking. However, 
for Fit la, although the HD 82943 system can be stable at least for 10 7 yr, it is still difficult 
to secure the resonant capture with N-body models due to the larger eccentricities of 0.41- 
0.54 in Kley et al. (2004), the origin for such configuration may attribute to an additional 
passing star's influence or pure gravitational scattering mechanism (Lin & Ida 1997; Adams 
& Laughlin 2003), to pump up the eccentricities. One of the final fates in the scattering 
procedure is that the inner planet is left alone in the system with an eccentric orbit while 
the outer planet is ejected (Ford et al. 2001) into the stellar space, nearly immune to the 
gravitational influence from the parent star. If the scattering does happen in the course of 
the migration, a casual rapid capture may be required to prevent one of the planets from 
throwing away from the system and subsequently two orbits are well locked into 2:1 resonance 
(see §3.2.1) at larger eccentricities. Alternatively, another possible scenario is that both of 
the eccentricities are excited to moderate values when the 2:1 resonance is being crossed 
(Chiang et al. 2002), subsequently the continuous resonance crossings may terminate for the 
suitable eccentricities evolved into the resonance geometry. 



6. Summary and Remarks 

In this paper, we have explored the stable geometry for a system with two planets 
involved in a 2:1 MMR, and we mainly concentrate on the study of the HD 82943 system by 
adopting two different sets of the orbital parameters. In our numerical simulations, we found 
there are three possible stable configurations for HD 82943 system: (l)Type I, only 9\ pa 0°, 
(2)Type II, 9 1 « 6 2 « 9 3 « 0° (aligned case), and (3)Type III, 6 1 « 180°, 9 2 « 0°, 6 3 w 180° 
(antialigned case). The direct integrations show that all stable orbits are related to the 2:1 
commensurability that remains the semi-major axes for two companions slightly vary about 
the initial values over the secular dynamical evolution for 10 7 yr, further the apsidal phase- 
locking for two planets can enhance the stability for this system, because the eccentricities are 
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simultaneously preserved to prevent the planets from frequent close encounters. Additionally, 
we comment that other systems trapped in 2:1 MMR (e.g. GJ 876 or possibly HD 160691) 
can select one of the aforementioned three stable topology in their authentic orbital motions, 
which are determined by their initial parameters. And we underline that such configurations, 
making up the dynamical families, can serve as a general regulation for the long-term stable 
orbits of the 2:1 resonant exosystems. In final, we summarize some conclusions: 

Using the earlier best-solutions of Fit 1, we found three types of stable orbits for HD 
82943, while we observed the steady configurations of Types I and II by employing the new 
solution of Fit 2 supplied by Mayor et al. (2004). In the meantime, we also proposed a 
semi-analytical model to study the — Aro Hamiltonian contour determined by the ini- 
tial parameters, and the closest level curves encompass the direct numerical results, which 
presents a good agreement between them. The theoretical figures can still provide valuable 
information on the dynamics for two planets. 

Given the stable geometry based on Fit 2, we then extensively examined the dependence 
of the stability of HD 82943 in the orbital parameter space and the planetary mass ratios. In 
the case of the non-coplanar circumstances, we found that stability requires that the relative 
inclination be ~ 25° or less. For a fixed planetary mass ratio ~ 1, the stable orbits for 
HD 82943 requires sini > 0.50. For the non-constant mass ratio case (where m-i is kept), 
the requirement of the stability is m\jm2 < 2, which indicates an upper boundary for the 
mass of the inner planet. Concerning the eccentricities, the system can be always steady 
when < e2 < 0.24 and < e\ < 0.60. In a word, these outcomes do demonstrate that 
the stability for HD 82943 is strongly sensitive to the sound planetary masses, rather lower 
relative inclination between two orbits, the eccentricities and other orbital parameters. In 
addition, we showed that the assumed terrestrial bodies cannot exist in the habitable zones 
of HD 82943 thanks to strong perturbations induced by two massive planets, but the Earth- 
like planets can be dynamically habitable in the GJ 876 system at ~ 1 AU in the numerical 
surveys. 

However, as the additional measurements will improve the present best-fit solution for 
HD 82943 with a corrected model by taking into account the planet-planet interactions. 
Thus, we hope to have more precise orbital solution not only to check for our presented 
predictions in this paper but to better make sense of the complete dynamics and the origin 
of this system in the near future. 

We thank G. Laughlin, M.H. Lee, S.J. Peale, J.C.B. Papaloizou and E. Thommes for 
informative discussions and helpful explanations, K. Gozdziewski for good comments. We are 
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A. The form of the total angular momentum with use of Jacobi coordinates 

The Hamiltonian with use of Jacobi coordinates has the following form: 

1 . 2 1 . 2 Gm Q mi Gm m 2 Gm x m 2 , . , , 
F = -<7iV + -a 2 r 2 2 , Al) 

2 2 r i r 02 r i2 

where 

m m 1 m 2 (m + mi) 

(7! = ■ , d 2 = ■ ■ . (A2) 

m + mi mo + mi + m 2 

Now we have two choices for the unperturbed Hamiltonian: 

A) 

F = 1<W - + 1^2 _ (A3) 

2 n 2 r 2 

1 2 Gmomi 1 2 Gm 2 (m +wi) 

^0 — ~<?l r l r -<7 2 f2 • l A4 J 

2 ri 2 r 2 

The difference between A) and B) is the last term. The choice A) is recently popular. In 
the treatment of the triple star system the choice B) is always adopted, since the form of 
the disturbing function becomes simpler for the case r\ < r 2 . 



B) 



A) 
B) 



The unperturbed Hamiltonian takes the following form for each case: 

Gm mi Gm m 2 , . 

F ° = —2^ 2^' (A5) 

_ Gm mi _ G(m + mi)m 2 ,^ , 

°~ 2ai 2a 2 1 ) 

The form of the total angular momentum J is 



A) 



B) 



[77< \ 7-, 9^ fTy m Q + mi+ m 2 ~ ~ 

cr 1 \ G(m + mi)ai(l - ef) + a 2 \ Gm ■ a 2 (l - e A 2 ), (A7) 

v V m o + m i 

ai^G(m + mi)ai(l - ef) + cr 2 \J G(m + m x + m 2 )a 2 (l - e|). (A8) 

In the above expressions the meaning of the osculating elements a 2 , e 2 is different between A) 
and B), because the unperturbed form is different. However, for the averaged Hamiltonian 
after the elimination of the short periodic terms, the difference due to the choice between 
the Jacobi coordinates and heliocentric (or astrocentric) coordinates is the order of the 
perturbation (mi /mo, m 2 /m ), which could be negligible. 
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Table I. The orbital parameters of HD 82943 planetary system (Fit l a ). m = 1.05M ( 





Fit 


la 


Fit 


lb 


Parameter 


Outer (b) 


Inner (c) 


Outer (b) 


Inner (c) 


m smi(M Jup ) 


1.63 


0.88 


1.63 


0.88 


P(days) 


444.6 


221.6 


444.6 


221.6 


a(AU) 


1.16 


0.73 


1.16 


0.73 


e 


0.41 


0.54 


0.41 


0.54 


fi(deg) 


14.79 


154.23 


327.65 


239.11 


cu(deg) 


96.0 


138.0 


96.0 


138.0 


Mean Anomaly(deg) 


97.86 


10.02 


53.40 


188.50 



a The orbital parameters for the planetary masses, orbital periods, semi- 
major axes, eccentricities and apsidal arguments are taken from the Geneva 
website, as of July 31, 2002. Both of the inclinations are assumed to be 0.5°. 
Note Fit la and Fit lb are, respective, one of Types III and II, amongst 
100 first runs. 
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Table 2. The orbital parameters of HD 82943 planetary system (Fit 2 a ). m = 1.15M ( 



Parameter 


Outer (b) 


Inner (c) 


m smi(M Jup ) 


1.84 


1.85 


P(days) 


435.1 


219.4 


a(AU) 


1.18 


0.75 


e 


0.18 


0.38 


fi(deg) 


120.91 


315.60 


cu(deg) 


237.0 


124.0 


Mean Anomaly(deg) 


168.56 


185.12 



a The orbital parameters for the planetary 
masses, orbital periods, semi- major axes, ec- 
centricities and the periastron arguments are 
adopted from Mayor et al. (2004), as of Oc- 
tober 13, 2003. The inclinations are the same 
as given in Table 1. The mean anomalies are 
derived at the epoch JD 2,452,396.82. Fit 2 is 
one of Type II of 100 second runs. 
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111 Known Companions, Aug 5, 2003 
http://exoplanets.org/almanacframe.html 



LLL 



0.0 0.1 0.2 



0.3 0.4 0.5 0.6 0.7 0.8 0.9 
Eccentricity 
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Fig. 1. — Left panel: the histogram for the eccentricities for 111 extrasolar planets, as of 
Aug. 5, 2003. Notice that over 50% of the planets have the eccentricities larger than 0.30, 
and HD 80606 b can occupy the eccentricity up to 0.93. Right panel: distribution of the 
semi-major axes plotted against the eccentricities. Note that 55 Cnc d can be distant as far 
as nearly 6 AU (see also Marcy et al. 2003) from its parent star. 
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Fig. 2. — Orbital evolution for Fit la (antialigned orbits). Upper panel: 9i librates (by 
dots) about 180° with a moderate amplitude of ~ 30°, 6 2 (by dots) about 0° with a small 
amplitude of ~ 10°, and 8 3 (by thick line) about 180° with a relatively small amplitude of 
~ 30° for t — 10 7 yr (to see more clearly, we simply display a snapshot for t = 2000 yr). 
Lower panel: a± and a 2 slightly vibrate about 0.73 and 1.15 AU for t = 10 7 yr. Notice that 
e\ and e 2 separately reside in (0.36, 0.56) and (0.39, 0.49). The 2:1 resonance is confirmed 
by the modulations of the semi-major axes and 6i (i = 1,2). Note that 9i and 6> 3 share 
the same libration period of ~ 300 yr. 
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Fig. 3. — Orbital evolution for Fit lb (aligned orbits). Upper panel: 9\ (by square sign) 
librates about 0° with a small amplitude of ~ 20°, but 9 2 (by plus sign) and 9 3 (by thick 
line) individually librate about 0° (a snapshot for t = 3000 yr ) with a large amplitude of 
~ 70°. In this resonant aligned geometry, we again observe that 9i, 9 2 and 6* 3 share the 
common librating periods of ~ 700 yr, coupled with those of the eccentricities of e\ and e 2 . 
Lower panel: a x and a 2 do not change dramatically but undergo small oscillations about 
0.73 AU and 1.15 AU for t = 10 7 yr, further the eccentricities e\ and e 2 librate in the range 
of (0.5, 0.8) and (0, 0.45), respectively, in the behavior of the converse cycles, owing to the 
conservation of the total angular momenta for the system. 



-39- 




Fig. 4. — e-Aw Hamiltonian contour map for Fit 1. Left panel: (a)For Fit la, the contour 
levels are exhibited by thin line and the numerical solutions shown by thick line (the inner- 
most curve about plus). The marked plus denotes the stationary solutions - (180°, 0.45) 
and (180°, 0.46), respectively. Right panel: (b)For Fit lb, thin line for contour levels and 
dotted points for numerical solutions, the librating centers are, respectively, (0°, 0.68) and 
(0°, 0.29). The e { — Aw figures both show a good agreement between the numerical results 
and the semi-analytical solutions. 
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Fig. 5. — Orbital variations for Fit 2 (aligned orbits). Upper panel: Q\ (by square sign) 
librates about 0° with a moderate amplitude of ~ 45°, but 6 2 (by plus sign) and 9 3 (by 
thick line) individually librate about 0° (a snapshot for t = 3000 yr ) with large amplitude 
of ~ 80°. Lower panel: the semi-major axes a\ and a 2 are almost unchanged and they 
modulate about 0.75 AU and 1.18 AU with relatively smaller amplitude for t = 10 7 yr, at 
this time, the amplitude of the oscillations for e\ and e 2 are not so large and they are just 
wandering in the span (0.34, 0.44) and (0, 0.25), respectively. Notice that the eccentricities 
have the libration periods of ~ 600 yr, which are coupled with those of 9\ and 9 3 . 




Fig. 6. — e-Azu Hamiltonian contour map for Fit 2. Same as the definition in Figure 4b. 
The equilibrium points are: (0°, 0.40) and (0°, 0.13), respectively. 
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Fig. 7a- 7b. — (a-b)The stable maps in the parameter spaces of [01,02] and [ei,e2]. (c-d)The 
stable maps in the parameter spaces of [M!,M 2 ] and [u^u^]. The sign crosses indicate the 
unstable orbits and filled circles for stable orbits lasting for 10 6 yr. The intersection between 
the vertical and horizontal lines represents for the present solution from Fit 2. 
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Fig. 8. — Initial orbital parameters for 100 seed Earth-like planets. Upper panel: the dis- 
tribution of [e, a], where a G [0.96AU, 1.05AU] with the incremental step of 0.01 AU, and 
e G [0,0.1], indicating the initial near-circular orbits for testing the dynamical habitabil- 
ity. Lower panel: eccentricity versus inclination, where % G [0°,5°], showing that these seed 
planets are almost in the same plane with the resonant pairs. 
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Fig. 9. — Left panel: (a)The numbers of the ejected terrestrial planet vs ejecting timescale. 
Notice that 95% of the orbits are ejected at r < 1.0 x 10 5 yr. Right panel: (b)A typical 
orbital evolution for the ejected orbits in the HD 82943 system: the semi-major axe a grows 
from ~ 1 AU to over 100 AU, the eccentricity e undergoes a rapid increase from ~ 0.1 to 1 
and the inclination i is also excited to a high value of ~ 40°. Finally, the assumed Earth-like 
planet is ejected away from the two-planet system at 4, 800 yr. 
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Fig. 10. — Orbital variations for possible Earth-like planet in the GJ 876 system. Both the 
semi-major axis a and the eccentricity e execute small fluctuations about 1 AU and 0.06, 
respectively, and the inclination % also remains less than 2 degrees over the same time span. 
And there is no sign to indicate that such regular orbits at ~ 1 AU with lower eccentricities 
will become chaotic for much longer time, even for the age of the star. 



